#version 120

#define    MULTIPLICATION          0
#define    EXPONENTIAL             1
#define    DIVISION                2
#define    SUBTRACTION             3
#define    ADDITION                4
#define    POTENCY                 5
#define    SIGNALCHANGE            6
#define    COSINE                  7
#define    LN                      8
#define    X                       9
#define    Y                      10
#define    CONSTANT               11
#define    SINE                   12
#define    TANGENT                13

#define    MAXARRAY              100

uniform int   expression[MAXARRAY];
uniform float constants[MAXARRAY];

//int expression[3] = int[](MULTIPLICATION, X, X);
//float constants[1] = float[](0.0);

uniform float xRange;
uniform float yRange;

struct AutoDiff
{
  float _   ;
  float dx  ;
  float dy  ;
  float dxx ;
  float dxy ;
  float dyy ;
};
AutoDiff new_autodiff()
{
  AutoDiff f ;
  f._  = 0.0 ;
  f.dx = 0.0 ;
  f.dy = 0.0 ;
  f.dxx = 0.0;
  f.dxy = 0.0;
  f.dyy = 0.0;
  return f ;	
}

AutoDiff new_autodiff_cste(float a)
{
  AutoDiff f ;
  f._  = a ;
  f.dx = 0.0 ;
  f.dy = 0.0 ;
  f.dxx = 0.0;
  f.dxy = 0.0;
  f.dyy = 0.0;
  return f ;
}
AutoDiff new_autodiff_varx (float x)
{
  AutoDiff f ;
  f._  = x ;
  f.dx = 1.0 ;
  f.dy = 0.0 ;
  f.dxx = 0.0;
  f.dxy = 0.0;
  f.dyy = 0.0;
  return f ;
}
AutoDiff new_autodiff_vary (float y)
{
  AutoDiff f;
  f._  = y ;
  f.dy = 1.0 ;
  f.dx = 0.0 ;
  f.dxx = 0.0;
  f.dxy = 0.0;
  f.dyy = 0.0;
  return f ;
}
AutoDiff ADsig(AutoDiff f)
{
  AutoDiff g ;
  g._ = -f._ ;
  g.dx = -f.dx ;
  g.dy = -f.dy ;
  g.dxx = -f.dxx;
  g.dxy = -f.dxy;
  g.dyy = -f.dyy;
  return g;
}

AutoDiff ADadd(AutoDiff f, AutoDiff g)
{
  AutoDiff z;
  z._   = f._   + g._  ;
  z.dx  = f.dx  + g.dx ;
  z.dy  = f.dy  + g.dy ;
  z.dxx = f.dxx + g.dxx;
  z.dxy = f.dxy + g.dxy;
  z.dyy = f.dyy + g.dyy;
  return z ;
}
AutoDiff ADsub(AutoDiff f, AutoDiff g)
{
  AutoDiff z;
  z._   = f._   - g._  ;
  z.dx  = f.dx  - g.dx ;
  z.dy  = f.dy  - g.dy ;
  z.dxx = f.dxx - g.dxx;
  z.dxy = f.dxy - g.dxy;
  z.dyy = f.dyy - g.dyy;
  return z ;
}
AutoDiff ADmul(AutoDiff f, AutoDiff g)
{
  AutoDiff z ;
  z.dxx = f._ * g.dxx + 2.0*(f.dx*g.dx) + f.dxx * g._;
  z.dxy = f._ * g.dxy + f.dx*g.dy + f.dy*g.dx + f.dxy * g._;
  z.dyy = f._ * g.dyy + 2.0*(f.dy*g.dy) + f.dyy * g._;
  z.dx  = f._ * g.dx  + f.dx * g._ ;
  z.dy  = f._ * g.dy  + f.dy * g._ ;
  z._   = f._ * g._ ;
  return z;
}
AutoDiff ADdiv(AutoDiff f, AutoDiff g)
{
  AutoDiff z ;
  z.dxx = (f.dxx*g._*g._ - f._*g.dxx*g._ -2.0*f.dx*g._*g.dx - 2.0*f._*g.dx*g.dx)/(g._*g._*g._);
  z.dxy = (f.dxy*g._*g._ - f.dy*g.dx*g._ - f._*g.dxy*g._ - f.dx*g._*g.dy + 2.0*f._*g.dx*g.dy)/(g._*g._*g._);
  z.dyy = (f.dyy*g._*g._ - f._*g.dyy*g._ -2.0*f.dy*g._*g.dy - 2.0*f._*g.dy*g.dy)/(g._*g._*g._);
  z.dx  = f.dx * g._  - f._ * g.dx ;
  z.dy  = f.dy * g._  - f._ * g.dy ;
  z._   = f._ / g._ ;
  return z;
}

AutoDiff ADpow(AutoDiff f, float r )
{
  AutoDiff z ;
  if(r == 0.0)
  {
    z._   = 1.0;
    z.dx  = 0.0;
    z.dy  = 0.0;
    z.dxx = 0.0;
    z.dxy = 0.0;
    z.dyy = 0.0;
  }
  else
  {
    z._   = pow(f._,r) ;
    z.dx  = r*pow(f._,r-1.0)*f.dx ;
    z.dy  = r*pow(f._,r-1.0)*f.dy ;
    z.dxx = r*( (r-1.0)*pow(f._,r-2.0)*f.dx*f.dx + pow(f._,r-1.0)*f.dxx );
    z.dxy = r*( (r-1.0)*pow(f._,r-2.0)*f.dx*f.dy + pow(f._,r-1.0)*f.dxy );
    z.dyy = r*( (r-1.0)*pow(f._,r-2.0)*f.dy*f.dy + pow(f._,r-1.0)*f.dyy );
  }
  return z ;
}

AutoDiff ADexp(AutoDiff f)
{
  AutoDiff z ;
  z._   = exp(f._) ;
  z.dx  = f.dx * z._ ;
  z.dy  = f.dy * z._ ;
  z.dxx = f.dxx * z._ + f.dx * z.dx ;
  z.dxy = f.dxy * z._ + f.dy * z.dx ;
  z.dyy = f.dyy * z._ + f.dy * z.dy ;
  return z ;
}
AutoDiff ADlog(AutoDiff f)
{
  AutoDiff z ;
  z._   = log (f._) ;
  z.dx  = f.dx /f._ ;
  z.dy  = f.dy /f._ ;
  z.dxx = f.dxx /f._ - z.dx * z.dx ;
  z.dxy = f.dxy /f._ - z.dx * z.dy ;
  z.dyy = f.dyy /f._ - z.dy * z.dy ;
  return z ;
}

AutoDiff ADcos(AutoDiff f)
{
  AutoDiff z ;
  z._  = cos(f._) ;
  z.dx  = -f.dx * sin(f._) ;
  z.dy  = -f.dy * sin(f._) ;
  z.dxx = -f.dxx * sin(f._) - f.dx * f.dx * z._  ;
  z.dxy = -f.dxy * sin(f._) - f.dx * f.dy * z._  ;
  z.dyy = -f.dyy * sin(f._) - f.dy * f.dy * z._  ;
  return z ;
}
AutoDiff ADsin(AutoDiff f)
{
  AutoDiff z ;
  z._  = sin(f._) ;
  z.dx  = f.dx * cos(f._) ;
  z.dy  = f.dy * cos(f._) ;
  z.dxx = f.dxx * cos(f._) - f.dx * f.dx * z._  ;
  z.dxy = f.dxy * cos(f._) - f.dx * f.dy * z._  ;
  z.dyy = f.dyy * cos(f._) - f.dy * f.dy * z._  ;
  return z ;
}
AutoDiff ADtan(AutoDiff f)
{
  AutoDiff z ;
  z._   = tan(f._) ;
  z.dx  = f.dx * (1.0 + z._*z._) ;
  z.dy  = f.dy * (1.0 + z._*z._) ;
  z.dxx = f.dxx * (1.0 + z._*z._)  + f.dx * 2.0 * z.dx * z._ ;
  z.dxy = f.dxy * (1.0 + z._*z._)  + f.dx * 2.0 * z.dy * z._ ;
  z.dyy = f.dyy * (1.0 + z._*z._)  + f.dy * 2.0 * z.dy * z._ ;
  return z ;
}

int numberOfOperands(in int o)
{
    if(     o == EXPONENTIAL    || o == SIGNALCHANGE || o == COSINE      || o == LN       || o == SINE   || o == TANGENT)
        return 1;
    if(o == MULTIPLICATION || o == DIVISION     || o == SUBTRACTION || o == ADDITION || o == POTENCY)
        return 2;

    return 0;
}

AutoDiff opera(int o, AutoDiff first, AutoDiff second)
{

    if( o == EXPONENTIAL)
        return ADexp(first);
    if( o == SIGNALCHANGE)
        return ADsig(first);
    if( o == COSINE)
        return ADcos(first);
    if( o == SINE)
        return ADsin(first);
    if( o == LN)
        return ADlog(first);
    if( o == TANGENT)
        return ADtan(first);

    if( o == MULTIPLICATION)
        return ADmul(first , second);
    if( o == DIVISION)
        return ADdiv(first , second);
    if( o == SUBTRACTION)
        return ADsub(first , second);
    if( o == ADDITION)
        return ADadd(first , second);
    if( o == POTENCY)
        return ADpow(first, second._);

    return first;
}

struct PILHA_O
{
    int o[MAXARRAY];
    bool hasFirst[MAXARRAY];
    int topo;
};

struct PILHA_R
{
    AutoDiff f[MAXARRAY];
    int topo;
};

PILHA_O po;
PILHA_R pr;

void pusho(int o)
{
    po.o[po.topo] = o;
    po.hasFirst[po.topo++] = false;
}
void rePusho(int o)
{
    po.o[po.topo] = o;
    po.hasFirst[po.topo++] = true;
}
void popo(out int o, out bool hf)
{
    o = po.o[--po.topo];
    hf = po.hasFirst[po.topo];
}
AutoDiff evaluate(AutoDiff x, AutoDiff y)
{
    po.topo = 0;
    pr.topo = 0;
    int idx = 0;
    int constIdx = 0;
    bool getNext = true;
    AutoDiff f1 = new_autodiff();
    AutoDiff f2 = new_autodiff();
    bool hf;
    int atual;

    do
    {
        if(!getNext)
        {
            popo(atual,hf);
            if(numberOfOperands(atual) > 1 && !hf)
            {
                rePusho(atual);
                getNext = true;
            }else
            {
                if(numberOfOperands(atual) > 1)
                {
                    f2 = pr.f[--pr.topo];
                    f1 = pr.f[--pr.topo];
                }
                else
                {
                    f1 = pr.f[--pr.topo];
                }
                pr.f[pr.topo++] = opera(atual,f1,f2);
            }
        }else
        {

            atual = expression[idx++];

            if(atual == CONSTANT)
            {
                pr.f[pr.topo++] = new_autodiff_cste(constants[constIdx++]);
                getNext = false;
            }
            else if(atual == X)
            {
                pr.f[pr.topo++] = x;
                getNext = false;
            }
            else if(atual == Y)
            {
                pr.f[pr.topo++] = y;
                getNext = false;
            }
            else
            {
                pusho(atual);
            }
        }
    }while(po.topo != 0);
    return pr.f[--pr.topo];
}

AutoDiff f(float x, float y)
{
  return evaluate(new_autodiff_varx(x),new_autodiff_vary(y));
}
float curvature( AutoDiff f )
{
  float denom = sqrt( f.dx*f.dx + f.dy*f.dy ) ;

  if(abs(denom) < 0.00001)
      return 0.0;

  return
  - (f.dxx * f.dx * f.dx
     - 2.0 * f.dxy * f.dx * f.dy
     + f.dyy * f.dy * f.dy ) /
  ( denom * denom * denom ) ;
}

void shading(vec3 n, vec3 p, inout vec3 color)
{
    vec4 lp = gl_LightSource[0].position;
    vec3 l = normalize(lp.xyz - p);
    float diff = max(0.0, dot(n, l));
    float spec = 0.0;
    if(diff > 0.0)
    {
        vec3 v = normalize(-p);
        vec3 r = normalize(reflect(-l,n));
        spec = pow(max(0.0, dot(v,r)) , 64.0 );
    }

    color = color*0.5 + color*diff*0.5 +  spec;
}

void main()
{
	vec4 v = gl_Vertex;
	v.x *= xRange;
	v.y *= yRange;
	AutoDiff resp = f(v.x, v.y);
	v.z = resp._;

        vec3 color = vec3(curvature(resp)*0.05, 0.3, 0.3);

        vec3 n = gl_NormalMatrix * vec3(-resp.dx, -resp.dy, 1.0);
        n = normalize(n);

        vec4 eyeCoord = gl_ModelViewMatrix * v;

        shading(n, eyeCoord.xyz , color);

        gl_Position = gl_ProjectionMatrix * eyeCoord;
        gl_FrontColor = vec4(color, 1.0);
}
